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New capabilities have been developed for a Navier-Stokes solver to perform steady-state simulations more effi- 
ciently. The flow solver for solving the Navier-Stokes equations is based on a combination of the lower-upper fac- 
tored symmetric Gauss-Seidel implicit method and the modified Harten-Lax-van Leer-Einfeldt upwind scheme. 
A numerically stable and efficient pseudo-time-marching method is also developed for computing steady flows over 
flexible wings. Results are demonstrated for transonic flows over rigid and flexible wings. 


Introduction 

T HIS paper summarizes new developments in the static option of 
the NASA Ames Research Center’s aeroelasticity simulation 
code ENSAERO, which computes steady-state flowfields and static 
deformations of structures simultaneously. This code is capable of 
computing aeroelastic responses by simultaneously integrating the 
Euler/Navier-Stokes equations and the modal structural equations 
of motion using aeroelastically adaptive dynamic grids. 1 It was en- 
hanced with the upwind option and applied to transonic flows from 
small to moderately large angles of attack for lighter wings undergo- 
ing unsteady motions. 2 - 3 Next, it was extended to simulate unsteady 
flows over a wing with an oscillating trailing-edge flap 4 The geo- 
metric capability of the code has further been extended to handle 
a full-span wing-body configuration with control surfaces. 5 As an 
option, the structure can be modeled using shell/plate finite element 
formulation. 6 

The present research has been performed to improve the static 
aeroelastic option of the code. The computational fluid dynamics 
part of the code has been completely rewritten with new algorithms. 
These new techniques do not change the accuracy of the code, but 
they make the code more efficient and robust. The lower-upper 
factored symmetric Gauss-Seidel (LU-SGS) method 7 is employed 
to reduce the arithmetic operation count in the implicit solver 
The modified Harten-Lax-van Leer-Einfeldt (HLLE) upwind 
scheme 8 is used to obtain a robust flow solver with reasonable costs. 
The code originally employed the streamwise upwind algorithm as 
an upwind option. This algorithm brings multidimensional informa- 
tion to the upwind technique and gives better accuracy than conven- 
tional dimensional-split upwind techniques. After the higher order 
extension is made in a standard manner, however, the improvements 
in accuracy become small, while the additional computational cost is 
high. Although it is acceptable for a research code to obtain the best 
accuracy on a given grid, efficiency is also important for a produc- 
tion code. Therefore, the streamwise upwind algorithm is replaced 
with the dimensional-split upwind algorithm. In addition, the HLLE 
upwind scheme is a positively conservative scheme. 9 The present 
modified HLLE scheme has the same resolution, but it is more robust 
than the widely used Roe upwind scheme. 10 Its arithmetic operation 
count still remains comparable to the Roe scheme. 
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In aeroelastic calculations, structural oscillations are often ob- 
served during the initial transient period, which leads to divergence 
of flow calculations or slow convergence. To reduce such oscilla- 
tions, the structural dynamics part of the code is also modified to 
use artificial structural damping along with a pseudo-time-marching 
method. 

To demonstrate the capability of the resulting code, three test cases 
are shown: computation of a flat-plate boundary layer, transonic 
flow around the ONERA M6 wing, and transonic flow about an 
aeroelastic wing. 6 11 


Algorithm Development 

ENSAERO computes the three-dimensional Euler or thin-layer 
Navier-Stokes equations with the Baldwin-Lomax model for the 
fluid flow part. The description of those equations and the turbulence 
model can be found in Refs. 1 and 12, for example. The essence ot 
the new algorithms can be described in the one-dimensional equa- 
tions as shown in the following sections. 

Modified HLLE Scheme 

The conservation form of the one-dimensional Euler equations is 
Qi + F x =0 (la) 

where the conserved quantities Q and flux F are 


/ p \ 

PU 

P« , F = 

pu 2 + p 

W 

_ (e + p)u __ 


and where p is the density, u is the velocity, and e is the total energy 
per unit volume. The pressure p is related to the conserved quantities 
through the equation of state for a perfect gas 

P = (y - 1)(* - pu 2 / 2) (lc) 


The cell interface flux F LR can be evaluated by the HLLE scheme 9 
as 


where 


F LR = {(F L + F R -RALAQ) 

(2a) 
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and where R , A, and L are the right eigenvector, eigenvalue, and 
left eigenvector matrices of the Roe-averaged Jacobian, 10 / is the 
identity matrix, and the subscripts L and R indicate the left and right 
states. (The overbar indicates Roe-averaged quantities.) The HLLE 
scheme defines b# and b~ L as 


at the cell interface. 14 For example, the interpolated pressure are 
given as 

Pu+\ — 1 1 + ^[(1 ~ -f (1 + k) A] 1/?, (6a) 


bx ~ max(w F c, w* F c R , 0) 

(3a) 

b~ — min(w — c t u L — c Lt 0) 

(3b) 


This scheme satisfies all of the stability, entropy, and posi- 
tively conservative conditions required for the nonlinear difference 
equations. 9 

The HLLE scheme approximates the solution of the Riemann 
problem with two waves propagating with speed of bn = max(u F 
c, u R F c R ) and b L — min(u — c, Ul — Cl) and a state Q LR between 
those waves. Compared with the Roe scheme, the HLLE scheme 
introduces large numerical dissipation at contact discontinuities. A 
newly proposed modification 8 improves the resolution at contact 
discontinuities by replacing A = diag[Aj , X 2 > A3] in Eq. (2a) with 

A = diag[A ( — 25min(6j, b l), \ 2 , A3] (4a) 


P/ti+l = { 1 - 7KI 4- k)V F (1 — k)A] \pi 


7 + t 


(6b) 


where V and A are backward and forward difference operators, 
respectively. For a one-dimensional or Cartesian grid case, third- 
order interpolation can be obtained from k = 1/3. One limiter 
function <p from Koren 15 is given as 

3Vp,--Api 

<P(Pi) = rr? ~ — (7) 

2(A p, - VpiY + 3 V pi ■ Api 

Since this limiter reacts to even small oscillations in smooth regions, 
convergence to steady state is often stalled. To avoid clipping smooth 
extrema, a modification to the limiter was proposed in Refs. 16 
and 1 7 as 


Hpd = 


3 V pi • Ap, + e.f 

2(A pi - Vp.Y + 3V/>, ■ Ap, + ej 


( 8 ) 


where 



and where cr\ = Ap — A p/c 2 . Thus 5 = 1/2. The resulting scheme 
reduces to the Roe scheme when 8 = \f2{pi R /\o\\ — ► ooas|<7i| — ► 
0), and to the scheme when 8 — 0(p LR /\cr\\ -> 0 as |ai| -> co). 
Because ctj represents a jump in entropy, it is zero for isentropic 
flows. Then the present scheme results in the Roe scheme. As the 
jump in entropy becomes large, the present scheme turns into the 
HLLE scheme. 

To compute the present modified HLLE flux, Eq. (2), there is no 
need to do a matrix computation as suggested for the Roe scheme. 13 
In the actual implementation, Eq. (2) is rewritten as 


F L r 


where 
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Again by replacing Aj with A) by using Eq. (4), the present algorithm 
is obtained. The extension to the three dimensions is straightforward 
by dimensional splitting. 


Modified Differentiable Limiter 

Higher order numerical fluxes are obtained from higher order 
interpolation of the primitive variables for the left and right states 


In this paper, the threshold sf is given by ef — max(3.0|Vf |. 3 , 
10~ 12 ) where £ denotes the generalized coordinate. 


LU-SGS Method 

Discretizing Eq. (1), the LU-SGS method 7 can be described as 
{/ + fc[Xp(A),/ - A+JK/ +h X p(A)I];' 

x {/ + h[ XP (A),l + Ar +1 ]}Afi = (9a) 

where h — A/ /Ax, A = dF/dQ, p(A) denotes the spectral radius 
of A, 

/C = \[A±xp(A)I)} (9b) 

and x = 1.01 typically. Its extension to three dimensions is straight- 
forward. Note that this is a two-factored scheme so that the algorithm 
can be written as 
Forward sweep: 

Ag; = [-/i(f; +i -F'^+hA'AQ’.'^/iX+hxpiAi)) 

(10a) 


Backward sweep: 

A Qi = Ag; - hA~ AQj+i/[\ F hxpiAi)) (10b) 


where 1 F hxp(Aj) is only a scalar quantity. Thus, it requires no 
block-matrix inversion. 

To vectorize these sweeps in three dimensions, a hyperplane is 
necessary. Suppose i, j t and k denote the indices of a grid point, a hy- 
perplane is expressed as i F j F& = const. The three-dimensional ar- 
ray (i, j, k) is then converted to two-dimensional array (/, m) where 
/ denotes the point in the hyperplane and m = i F j + k represents 
the hyperplane. This enjoys a very long vector length, but the per- 
mutation requires additional memory space. Instead, the existing 
three-dimensional array can be used as it is by vectorizing only in 
the i index. The outer DO loop is for m. The next loop is for k. The 
inner DO loop is for i y and j is determined as j = m — k — i. Its 
vector length is shorter than the former, but it requires no additional 
memory space for vectorization. 

The split Jacobian matrices A* mimic the upwind Jacobians. For 
a natural upwinding, A + = A and A~ = 0 when u > c. However, 
the definition of A* in Eq. (9b) gives A “ ^ 0. A simple modification 
to the split Jacobian was tried as 


A ± = y [A ± p(A)l] ± ( JLJl p(A) , 

where 

± _ I 1. if ± (« ± c) > 0 

5 _ jo, if±(«±c) <0 


(Ha) 


(Hb) 
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In the authors' experience, it gives slightly better convergence 
than the unmodified case. For the flat-plate boundary-layer case, 
however, it did not make the convergence any better or worse. 

For the thin-layer Navier-Stokes computations, the split 
Jacobians in the viscous direction are further modified as 

A ± = y [A ± P(A)1 ] ± ± yyy |Vf 1 2 / (12) 


where fx is the laminar viscosity and (X T the turbulent viscosity. 

In addition, Eq. (9b) introduces large dissipation through the time 
integration. Only when the solution converges to a steady state, its 
effect vanishes. However, the LU-SGS method which is first-order 
accurate in time similar to other implicit methods, can be used for 
unsteady calculations by using special methods. For example, the 
Newton iteration 1819 can be used for unsteady computations, which 
not only removes the excess dissipation but also gives the second- 
order accuracy in time. 

In the following calculations, a locally varying time step was 
taken as 


AtiJ.k — A ^global 


1 + 0.0005 


03) 


where J is the Jacobian of transformation and A/ g | oha i = 1-5, de- 
pending on the problem. The Jacobian scaling in the denominator 
has been used commonly (for example, see Ref. 20). The additional 
scaling in the numerator has been also used widely without doc- 
umentation. It prevents the time step from becoming too small as 
the Jacobian becomes very large. With the LU-SGS scheme, this 
Jacobian scaling for the time stepping might not be necessary since 
the scheme is unconditionally stable. However, it was used to make 
fair comparisons with the previous version of the code. 


Aeroelastic Equation of Motion 

The governing aeroelastic equations of motion of a flexible wing 
are solved using the Rayleigh-Ritz method. It is assumed that the 
deformed shape of the wing can be represented by a set of dis- 
crete displacements at selected nodes. From the modal analysis, the 
displacement vector {d} can be expressed as 

{<*} = [«(*} O 4 ) 

where [tp] is the modal matrix and {q} is the generalized displace- 
ment vector. The final matrix form of the aeroelastic equations of 
motion is 


{M}{q} + [D]{q) + [K){q\ = [Z) (15) 

where [A/], [£>], and [K] are modal mass, damping, and stiffness 
matrices, respectively. The term (Z) is the aerodynamic force vec- 
tor defined as pU^Y [L]{ AC r }/2, and [LI is the diagonal area 
matrix of the aerodynamic control points. The aeroelastic equation 
of motion, Eq. (15), is solved by a numerical integration technique 
based on the linear acceleration method. 21 

Since {q ( = {q ) = 0 at a steady state, a primary static option uses 

[K][q}= {Z} (16) 

However, this often causes significant numerical oscillations be- 
cause there is no damping in the system. These oscillations are 
initiated by numerical transients. To obtain a smooth numerical 
transition, a pseudo-time-marching method can be taken by using 
Eq. (15) where the time step is set as A/ stnjcUirc = 0.02Ar g | oba i typ- 
ically. The damping coefficients [D] are set so that all aeroelastic 
modes damp out quickly. The damping matrix is computed using 
the relation 


[D] = a[M] + fi[K] (16a) 

where a and fi are functions viscous damping coefficients of each 
mode. More details about the computation of [D] can be found in 
Chap. 9 of Ref. 22. In this work, viscous damping coefficient val- 
ues for all modes are set to a high value of 0.80. This large value 
for the damping coefficient gives an artificial “shock absorber" for 
structural oscillations caused by flow transients. 


Results 

Flat-Plate Boundary Layer 

First, performances of the new algorithms are examined. Figure 1 
shows a comparison of computed M-velocity profiles of a laminar 
layer on a flat plate. The freestream Mach number is 2, the Reynolds 
number is 1 x 10 6 , and the freestream temperature is 222 K. The 
Sutherland law was used to compute the laminar viscosity. This 
computation was performed on a 1 1 1 x 5 x 69 grid by enforcing 
two dimensionality through five planes. The grid was generated 
algebraically. The time step Ar g i oba i used was 2.6 for all cases 

The upwind (modified HLLE) solution coincides with the similar- 
ity solution perfectly, whereas the central difference (CD) solution 
gives slight smearing at the outer edge of the boundary layer. The 
CD option uses the combination of the second- and fourth-order 
dissipation with the pressure switch. 20 The smoothing coefficients 
were set to 0.25 and 0.005 for the second- and fourth-order dissipa- 
tion terms, respectively. The upwind algorithm performs better as 
expected. Since obtaining grid convergence is not practical for most 
cases, the upwind option is generally recommended. Figure 2 shows 
a comparison of the convergence history (in terms of the averaged 
residual for all five equations multiplied by the Jacobian) between 
unmodified and modified limiters. The modified limiter, Eq. (8), 
shows a clear improvement over the original limiter. 

Figure 3 shows a comparison of convergence history among four 
different schemes. The original ENSAERO has the CD option with 
the diagonal approximate factorized alternating direction implicit 
(ADI) solver 20 and the upwind option using the streamwise upwind 
algorithm with the LU-ADI solver. 2 The new version of the code 
employs the LU-SGS solver for both CD and upwind options. The 
computational time per grid point per iteration in microseconds was 
4.0 for the CD option with the LU-SGS solver, 6.0 for the CD option 
with the diagonal ADI solver, 9.5 for the modified HLLE upwind 
option with the LU-SGS solver, and 13.8 for the streamwise upwind 
option with the LU-ADI solver, on a single Cray C90 processor. It 
should be noted that there is a disadvantage in the vectorization for 
the ADI solver in this comparison because the grid has only five 
planes in one direction. 



Fig. X Comparison of computed laminar boundary-layer profiles on 
a flat plate between the central difference and the upwind algorithms, 
Afoo =2.0 and Re = 10*. 



Fig. 2 Comparison of convergence history between the original and 
modified differentiable limiters, = 2.0 and Re = 10 6 . 



OBAYASHI AND GURUSWAMY: STATIC AEROELASTIC COMPUTATIONS 


1137 


The LU-SGS solver shows better convergence for both CD and 
upwind options than the ADI solvers. The advantage of using the 
LU-SGS solver is more apparent when the residual is replotted 
against CPU time in seconds using a single Cray C90 processor as 
shown in Fig. 4. Because of fewer operations and better convergence 



Fig. 3 Comparison of convergence history among four methods using 
the present and the previous version of the code, M - 2.0 and Re = 10 . 



Fig. 4 Comparison of convergence history in CPU seconds among four 
methods using the present and the previous versions of the code, - 
2.0 and Re = 10*. 


rates, the LU-SGS solver is about four times faster than the diag- 
onal ADI solver with the CD option. Even the new upwind option 
using the modified HLLE scheme is faster than the original CD op- 
tion with the diagonal ADI solver. The original upwind option is 
obviously the slowest because of the multidimensional upwinding. 
To obtain further acceleration for the LU-SGS solver, one might 
consider implementing the multigrid technique. 23 


ONERA M6 Wing 

To further validate the flow solver, transonic flows over the 
ONERA M6 wing have been computed by using the new upwind 
option of the code. A series of experimental data can be found in 
Ref. 24. Following a recent report of computational results, 25 the 
two cases chosen were one for Mach number 0.8395 and angle of 
attack of 3.06 deg and the other for Mach number of 0.8447 and 
angle of attack of 5.06 deg. The Reynolds number was set to 1 1 
x 10 6 for both cases. A C-O grid of size 193 x 34 x 49 was used. 
This is the coarser grid used in Ref. 25. The nondimensional time 
step was set to 2.0. 

High angle of attack cases with a large flow separation present a 
severe test for turbulence models. Even the Baldwin-Lomax model 
was found to be sensitive for a particular function evaluation. 26 The 
model requires the definition of w max , (u ^ ]F = — “mm)’ w ^ich 

is supposed to be M max = u\ p^, not the maximum of u along the 
grid line normal to the wall. For the flat-plate boundary layer with 
a uniform outer flow, for example, it does not matter how u max is 
defined. However, u max should be defined at the F max location. 

The 5.06-deg angle of attack case of ONERA M6 wing was sen- 
sitive to this definition of w max . The computation didn’t converge 
by using the true definition of M max , whereas the computation did 
converge by using the nominal the maximum of u along the 
grid line normal to the wing surface. Since turbulent viscosity is 
increased by overestimating M maX) more numerical viscosity is in- 
troduced into the computation. It probably helped stabilize the com- 
putation. In the following calculations, the nominal w max is used in 
the Baldwin-Lomax model. 

Figure 5 shows the surface pressure distributions compared with 
the experiment at the 44, 65, 80, and 90% spanwise sections for 
the 3.06-deg case. The computational results agree reasonably well 
with the experiment except on the upper surface of the 80% spanwise 
section where a finer grid may be required to resolve the merger of 






x/c x/c 


Fig. 5 Comparison of computed pressure distributions on the ONERA M6 wing with experiment, = 0.84, Re c = 11 X 10*, and a. = 3 deg. 
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two shock waves. Figure 6 shows the comparison of the convergence 
history between the new upwind option (modified HLLE + LU-SGS) 
and the previous CD option (diagonal ADI) in terms of both the 
residual and lift coefficient. The new upwind option shows a better 
convergence history. It took roughly 4500 iterations to converge to 
three digits of accuracy in the lift coefficient for the new upwind 
option. TTie CD computation was performed with the time step 1.0, 
since it diverged with the time step 2.0. The smoothing coefficients 
in this case were also set to 0.25 and 0.0 1 for the second- and fourth- 
order dissipation terms, respectively, since the computation diverged 
with the fourth-order dissipation coefficient 0.005. The CD option 
required roughly 8700 iterations to converge to three digits in the 
lift coefficients. 

The computational time per grid point per iteration using the 
upwind option was 7.7 /xs on a C90 single processor. This number 
is better than the flat-plate case (9.5 /xs) because the grid is truly 
three dimensional. On the other hand, the diagonal ADI CD option 
requires 5.9 /xs. To obtain the three-digit convergence in the lift 
coefficient, the new upwind option required two-thirds of the total 
computational time of the previous CD option. Although the time 
per grid point per iteration is reasonably fast for both cases, the 
convergence is still slow because of the more complicated flow- 

Modified HLLE (LU-SGS) 



Fig. 6 Comparison of convergence history between the new upwind 
option and the previous CD option of the code, ONERA M6 wing, 

~ 0.84, Re c = 11 X 10 6 , and a = 3 deg. 


field in this case than that in the flat-plate case. To improve the 
convergence rate further, the multigrid technique may be useful as 
demonstrated on the same grid in Ref. 25. 

Figure 7 shows the corresponding surface pressure distribu- 
tions for the 5.06-deg case. The Baldwin-Lomax model predicts 
a stronger shock wave and less flow separation, similar to the re- 
sults in Ref. 25. 

Aeroelastic Wing 

To demonstrate static aeroelastic computations, a swept wing of 
aspect ratio 3 and taper ratio 1/7 with the NACA 65A006 airfoil 
section was selected. The planform of the wing is shown in Fig. 8. 
Experimental and numerical studies can be found in Ref. 1 1 and 6, 
respectively. A C-H grid of size 151 x 44 x 44 was used. The 
time step was set to 6.25 for all of the steady cases. Only the upwind 
(the modified HLLE and LU-SGS option) results are shown here. 
The computational time per grid point per iteration was 8.0 /xs for the 
fluid part and 0.7 /xs for the structure part, including computations 
for grid movement. 

In the following calculations, the freestream Mach number was 
set to 0.854, the Reynolds number was 0.597 x 10 6 , and the ratio of 
specific heats was 1.135 because the experimental fluid was Freon. 
The Baldwin-Lomax model was used to compute turbulent viscos- 
ity. The dynamic pressure was set to 0.7 psi. The natural vibra- 
tion modes of the wing were calculated by the finite element plate 
model. 6 The computed frequencies for first three modes are 21.8, 
78.1, and 126 Hz, and the corresponding measured values are 21.6, 
79.7, and 121 Hz. Following Ref. 22, the values obtained fora and 
p are equal to 137.5 and 2.912 x 10“ 3 , respectively. 

Figure 9 shows responses of the leading-edge displacement at the 
wing tip. The displacement is plotted in inches, and 5000 iterations 
correspond to about 0.05 s in the unsteady computations. Small gen- 
eralized accelerations were given as initial conditions for the struc- 
ture. Starting from the steady-state flow solution at 1-deg angle of 
attack, three aeroelastic computations were performed. Two of them 
are steady-state calculations using the locally varying time stepping 
in the fluid part. One of them uses the pseudo-time-marching method 
for the structure using Eq. (15) (denoted as “pseudodynamic”), and 
the other uses the static equation, Eq. (16) (denoted as “static”). 
The third case is an unsteady computation using the time-accurate 





Fig. 7 Comparison of computed pressure distributions on the ONERA M6 wing with experiment, - 0.84, Re c = 11 x 10 6 , and a = 5 deg. 
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Fig. 8 Planform of an experimental swept wing. 


option of the code where the time step is constant for both fluid 
and structure part. It uses high damping coefficients 6 in Eq. (15) 
(denoted as “dynamic with damping”). 

The two steady-state computations converged within 1000 iter- 
ations. The unsteady result shows that the solution approaches the 
same steady state slowly. The unsteady computation with the damp- 
ing was continued up to 10,000 iterations and confirmed to converge 
to a displacement of 0.279 in. However, this option is too slow to 
be used for a production code. Figure 10 shows the magnified view 
of Fig. 9 for the first 1000 iterations. Both of the steady-state cal- 
culations go through the initial transient in the first 200 iterations 
and converge to 0.279 (three digits) in about 1000 iterations. The 
main difference of the two is the behavior during the initial transient. 
The pseudodynamic case shows very smooth transient, whereas the 
static case shows significant oscillations. Although it damps out 
quickly in this case, such oscillations often cause numerical insta- 
bility in the fluid part. Thus, the pseudo-time-marching method is 
very favorable for static aeroelastic computations. 

Since the pseudo-time-marching method for the structure gives 
a smooth transient, and computations can be started from the 
freestream condition impulsively. Note that the previous static 
aeroelastic calculations required 1000 iterations in addition to com- 
puting the initial steady-state flow solution without the structural 
dynamics. Figure 1 1 shows responses of the leading-edge displace- 
ment at the wing tip for this impulsive start case. The initial transient 
of the 1 -deg angle of attack case was damped out in 2000 iterations. It 
converged to 0.279 (three digits) in about 3400 iterations. This con- 
verged value coincides with that obtained in Fig. 10. The freestream 
condition was also set to 3-deg angle of attack. The convergence 
was slower than that for the 1-deg case, but it reached steady state 
roughly at 3500 iterations. 

Figure 12 shows the corresponding responses of the lift coef- 
ficients. The lift coefficients behave similarly to the leading-edge 
displacements shown in Fig. 11. The rigid wing cases without any 
structural dynamics are included for the comparison. The compar- 
ison between the flexible and rigid cases indicates that there is no 
penalty for having structural dynamics on convergence. In fact, the 
convergence history is almost identical for aeroelastic and rigid 
cases to a certain degree as shown in Fig. 13. Since the aeroelas- 
tic option does not slow down the convergence, it is more efficient 
and much simpler to perform static aeroelastic simulations directly 
starting from freestream conditions. 

The aeroelastic case seems to fall into the limit cycle earlier than 
the rigid case, although the 3-deg cases go into the limit cycle before 
they depart from each other. For the fluid part, the threshold in Eq. (8) 
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Fig. 9 Comparison of computed responses of leading-edge displace- 
ments at wing tip with steady and unsteady options of the code, Moo = 
0.85, Re c = 0.6 x 10 6 , and a = 1 deg. 



Fig. 10 Close-up view of computed responses of leading-edge displace- 
ments at wing tip. 


Pseudo dynamic, 1 deg 



Fig. 1 1 Computed responses of leading-edge displacements at wing tip 
with the pseudodynamic options of the code, M 0 0 = 0.85, Re c = 0.6 x 10 , 
and a = 1 and 3 deg. 


Pseudo dynamic, 1 deg 

Pseudo dynamic, 3 deg 



Fig. 12 Comparison of computed lift responses for flexible and rigid 
wing cases, M 0 0 - 0.85, Re c = 0.6 X 10 6 , and a = 1 and 3 deg. 
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Fig. 13 Comparison of convergence history for flexible and rigid wing 
cases, Moo — 0.85, Re c - 0.6 x 10 6 , and a - 1 and 3 deg. 





x/c 

Fig. 14 Comparison of computed surface pressures between flexible 
and rigid wing cases. Moo - 0.85, Re c = 0.6 x 10 6 , and a = 3 deg. 

can be tuned case by case to obtain better convergence. However, it 
is simply more practical to accept such limit cycles in convergence 
history for general applications. In such cases, one may choose 
either lift coefficient or leading-edge displacement as an index for 
convergence. 

In Fig. 12, the flexible cases consistently show lower lift at both 1 
and 3-deg angles of attack. Figure 14 shows comparison of surface 
pressure coefficients between the flexible and rigid wing cases at 
3-deg angle of attack. The plots confirm the lower lift for the flexi- 
ble wing. Figure 15 illustrates the corresponding tip displacement. 
The wing tip is bent up, and its leading edge is twisted down. The 
aeroelastic deformation reduces the effective angle of attack locally 



Fig. 15 Comparison of wing tip locations between flexible and rigid 
wing cases, Moo = 0,85, Re c = 0.6 x 10 6 , and a. - 3 deg (the vertical 
scale is exaggerated). 



Rigid 


Fig. 16 Comparison of wing shapes between flexible and rigid wing 
cases, Moo - 0.85, Re c = 0.6 X 10®, and a = 3 deg. 

near the tip, which results in the lower lift. Figure 16 shows the 
overall view of the wing deformation. 

Conclusion 

New capabilities to compute static aeroelasticity more efficiently 
have been added to ENSAERO code. The test case of computing 
a boundary layer on a flat plate demonstrates that the combina- 
tion of the LU-SGS implicit method and the modified HLLE up- 
wind scheme gives a good compromise of accuracy and efficiency 
for solving the Navier-Stokes equations. This option computes a 
steady-state solution faster than the original central difference option 
of the code. The capability of this new flow solver has been demon- 
strated successfully for computing transonic flows over wings. 

The existing procedure for static aeroelastic option requires a 
good initial guess, specifically a converged steady-state flow solu- 
tion. In this study, it is found that aeroelastic computations can be 
started from freestream conditions with the use of a pseudo-time- 
marching method for the structural equations of motion. The test 
case of computing transonic flows over a flexible wing indicates 
that static aeroelastic simulations can be performed at computa- 
tional effort similar to the cost of single-discipline steady-state flow 
simulations. 
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